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Abstract 

We constrain mass, lifetime and contribution of a very slowly decaying Ultra Heavy Dark 
Matter (UHDM) by simulating the cosmological evolution of its remnants. Most of interac- 
tions which participate in energy dissipation are included in the numerical solution of the 
Boltzmann equation. Cross-sections are calculated either analytically or by using PYTHIA 
Monte Carlo program. This paper describes in detail our simulation. To show the importance 
of the distribution of matter in constraining WIMPZILLA jl] characteristics, we consider 
two extreme cases: a homogeneous universe, and a local halo with uniform distribution. In 
a homogeneous universe, the decay of the UHDM with a mass ^ 10^'^GeU and a lifetime as 
short as a few times of the age of the Universe, can not explain the flux of the observed Ultra 
High Energy Cosmic Rays (UHECRs) even if the whole Dark Matter (DM) is composed of 
a decaying UHDM. If our simple halo model is reliable, in a uniform clump with an over- 
density of ^ 200 extended to lOOkpc, the lifetime can be ~ 10 — IOOtq, again assuming 
that DM is a decaying UHDM. We also compare our calculation with observed 7-rays at 
E - 10"ey by EGRET and CASA-MIA hmit at ~ IQ^^eV. They are compatible with a 
UHDM with relatively short lifetime. 

PACS codes: 95.35.-Kd, 98.70. Sa, 13.85.Tp, 98.70.Vc, 95.30.Cq, 04.25.Dm 



^Present Address: 03, impasse de la Grande Boucherie, F-67000, Strasbourg, France. 
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1 First Encounter 



Particle Physics today is a zoo with plenty of wild species very difficult to trace or capture. They 
are in general known under phyla WIMPs, LSP, Axions, Higgs, Heavy Neutrinos, etc. Most of 
these species are potential candidates of the Dark Matter. 

Recently a new phylum called WIMPZILLA |^ has been added to this zoo. The only common 
characteristic of the members of this family is their enormous mass, close to GUT scale lO^^GeV 
and their presumed long lifetime, much larger than the age of the Universe. 
Theoretical motivations for existence of these particles is not very new. Since early 90s, some 
compactification scenarios in string theory have predicted composite particles (e.g cryptons) 
with large symmetry groups and M > lO^^GeV. M-theory |Q] provides better candidates if 
compactification scale is much larger than Standard Model weak interaction scale . Messenger 
bosons in soft supersymmetry breaking models (see for review) also can have close to GUT 
scale masses. Being composite and decaying only through non-renormalizable interactions or 
having discrete gauge symmetry can make these ultra heavy particles meta-stable. Paramet- 
ric resonance or vacuum fiuctuation |2| at the end of inflation can produce a large amount 
of UHDM and unitarity constraint on their mass |^] can be overcome if they have never been 
thermalized [|[. 

If there was not other motivation for existence of an ultra heavy meta-stable particle, it was 
just one of many predictions of Particle Physics models waiting for detection. However, the 



discovery of Ultra High Energy Cosmic Rays (UHECRs) by large Air Shower detectors [10] [p 



1 12 1 (For a review of UHECRs detection and observed properties see [^]) is an observational 
motivation. The predicted GZK cutoff |14] in the spectrum of CRs at energies around ~ 10^^ 
due to interaction with CMBR and IR photons restricts the distance to sources to less than 
20 - 30Mpc. At present there are about 600 events with E ~ W^^eV, 50 with E > 4 x W^^eV 



and 14 events with E > lO'^^eV |T|] ||] including one with E ~ lO^^eF |Tl|]. The UHECR 
spectrum has a local minimum around E ~ 10^^ but it rises again at higher energies. 
Composition of the primary particles |jl^ [|l8| can be estimated from the shower and its muon 
content maximum position and their elongation rate in the atmosphere. In spite of uncertainties 
and dependence on hadrons interaction models at high energies |1£], all analyses of the data are 
compatible with a composition change from iron nuclei to proton at E > 10^^ eV IH |l|] [||]. 
Based on theoretical arguments however, some authors suggest that events with highest energies 
can be produced by heavy nuclei |^ (But see also [^] for maximum fly distance of Fe nuclei). 
Some of ideas about the origin of UHECRs has been reviewed in [^] and references therein. 
For the sake of completeness here we briefly review arguments for and against conventional and 
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exotic sources. 

Classical Candidates: In a recent review of conventional candidates of UHECR sources p^ ], 
R. Blandford rules out practically all of them. 

It is expected that charged particles are accelerated in the shock waves of AGNs, SNs, or in- 
falling gas in rich galaxy clusters. Their maximum energy is proportional to the magnetic field 
strength and the size of the accelerator. Using this simple estimation, remnants of old SNs are 
expected to accelerate protons up to ~ lO^^eV and iron to ~ lO^^eV [24|. A ~ lOOTeV 7-ray 
is expected from synchrotron radiation of these protons but nothing has been observed [p^ j. 
Shock front of AGNs relativistic jets can accelerate protons to ~ 10^^ eV pq]. Central black 



hole in Fannaroff- Riley galaxies |27] or the remnant of QSO [£8[ can produce UHECRs up to 
E ~ lO^^eV with a marginal maximum energy of 4 x lO'^^eV |^^. Adiabatic expansion [^0| and 
in situ interaction with radiation and matter fields however reduce the achievable maximum 
energy by orders of magnitude. Abrupt termination of the acceleration zone or composition 
change from p to n [30| increases the chance of particles to keep their energies. The former 
however needs a fine tuning of the source structure and in the latter case particles lose part of 
their energy anyway. 



Gamma Ray Bursts (GBR) also have been proposed as a candidate source of UHECRs [31]. A 



simulation |32| of cosmological distribution of sources with a power law flux of UHECRs shows 
however that the expected flux on Earth is much lower than observed value. 
In ||3^ it is argued that Poisson noise in the process of proton interaction with background 
photons leaves a non-interacting tail in the flux of UHECRs and increases the probability of 
detecting UHECRs from further distances. Optical depth of protons around GZK cutoff can be 
roughly estimated by Topt ~ aricmb- For a ^ OAbmb close to the resonance, the probability of 
non- interacting in a distance of 30Mpc is at most ~ 10^^. 

Correlation with Astronomical Objects: At present no serious astronomical candidate 
source has been observed. A correlation between Super Galactic Plane and UHECR events 



direction has been claimed [38|, but ruled out by other analyses It is plausible that the 



apparent clustering of events reported by AGASA collaboration |15| originates from caustics 
generated by the galactic magnetic held Q p^] pO[. 



Most of UHECRs burst models predict only one important nearby source [31| [33|. M87 in Virgo 
Cluster is practically the only conventional candidate that respects the distance constraint and 
is probably able to accelerate protons to ultra high energies. Recently, it has been shown that 
galactic wind and its induced magnetic field can deflect protons with E > 10^^ eV if the magnetic 
field is as high as 7/iG. In this case, most of UHECR events point to Virgo Cluster [^] . Even 
with such strong magnetic fleld, this model can correlate most energetic events with M87/Virgo 
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only if the primaries are He nuclei. 

There are also claims of correlation between the direction of UHECRs and pulsars |41] or 



QSOs |42|. The latter is reliable only if ultra high energy i/s produce protons by interacting 



with relic or a halo of v |42]. Recent works [44| however rule out a large part of the parameter 
space. 

Close to uniform distribution of UHECRs events has been concluded to be the evidence that 
UHECRs originate from some extra-galactic astronomical objects and not from a decaying 



UHDM in the Galactic Halo g|] |4|] g^]. MACHOs observation Q however shows that the 
Halo has a heterogeneous composition and a precise modeling of the anisotropy must take into 
account the distribution of various components. 

Exotic Sources: Ultra Heavy particles can be either long life DM or short life particles pro- 
duced by the decay of topological defects In the latter case, the heavy particles decay in 
their turn to ordinary species and make the observed UHECRs. The result of the decay of a 
heavy short life particle and a heavy relic can be quite the same, but their production rate and 
its cosmological evolution are very different and depend on the defect type and model [^ ]. 
The possibility that topological defects be the source of UHECRs has been studied exten- 



sively [pT| [m| |22]. Nonetheless, the observed power spectrum of LSS and CMB anisotropics 



rule out the existence of large amount of defects in the early Universe [^] [ 54 1 and consequently 
the chance that UHECRs be produced by their decay (there are also limits from high energy 
u 11). 

Another proposed source of UHECRs is the evaporation of primordial black holes (PBH). The 
Hawking temperature at the end of their life is enough high to produce extremely energetic ele- 
mentary particles like quarks and gluons and thus UHECRs. Most of models for production of 
PBH however needs fine tuning. Moreover, they are mainly produced when a large over density 



region crosses the horizon [55|. PBHs with present temperature of the same order as UHECRs 
energy must have an initial mass of 10^^ — W^^gr and thus formed when the temperature of the 
Universe was ~ lO^GeV. A thermal inflation [^ ] at EW scale would have reduced their density 
~ 10^^ times. In models with low reheating temperature, at these scales the Universe is not yet 
thermalized and parametric resonance and fluctuation production happens only at superhorizon 
scales p. 

Summarizing the discussion of this section, it seems that conventional sources of Cosmic Rays 
are not able to explain the observed rate of UHECRs j^^. Between exotic sources, the decay of 
a meta-stable ultra heavy particle seems to be the most plausible one. 
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2 WIMPZILLA Decay and Energy Dissipation of Remnants 



The decay of UHDM can have important imphcations for the evolution of high energy back- 
grounds. This can also be used for verifying this hypothesis and constraining the mass and 
lifetime of these particles. A number of authors have already tried to estimate the possible 
range of parameters as well as the flux of remnants on Earth (see |57] and [59| for defects, and 



1 58 1, |6C], 1 61] and |62| for UHDM). In the present work more interactions have been included 
in the energy dissipation of UHDM remnants. Hadronization is implemented using PYTHIA 
Monte Carlo program |^3| . We consider two distributions for UHDM. First we study the evolu- 
tion of the spectrum of stable particles i.e. e^,p^,i/, z/ and 7 in a homogeneous universe from 
photon decoupling to today. Then to show the effect of matter dumpiness, we simulate the 
Galactic Halo by simply considering a uniform over-density. A complete treatment of a clumpy 
universe will be reported elsewhere. We also estimate the effect of very slow decay of DM on 
the equation of state of the Universe and on the baryon and lepton asymmetries. 
In this section we describe the decay model of UHDM and interactions which are included in 
the simulation. 

2.1 Decay Model 

Theoretical predictions for mass and lifetime of UH particles cover a large range of values 
muH = 10^2 - W'^^eV and tuh = 10^ - lO^Oyr || g ||. Nevertheless, at the end of inflation 



it is more difficult to produce the highest range of the masses and special types of inflationary 
models [^] are needed. For UH particles make a substantial part of the Dark Matter today, 
their lifetime must be at least comparable to the age of the Universe. In this work we perform 
the simulation for mun = W'^'^eV and muH = W^'^eV. For lifetime, tjjh = 5ro and tuh = SOtq, 
where tq is the age of the Universe, are studied. These values are smaller than what have been 



used by other groups |58| |60| We show below that taking into account a realistic model for 
energy dissipation of remnants, even these relatively short lifetime can not explain the flux of 
UHECRs in a homogeneous universe. For some halo and IR background models, these lifetimes 
or slightly larger ones are compatible with observations. 

The decay modes of UHDM are very model dependent. It is very likely that they don't decay 
directly to known particles and their decay has a number of intermediate unstable states that 
decay in their turn. It is also very probable that remnants include stable WIMPs which are not 
easily observable. To study the maximal effects of the decay on high energy backgrounds, we 
assume that at the end, the whole decayed energy goes to stable visible particles. 



5 



Most of WIMPZILLA models consider them to be neutral bosons. Due to lack of precise in- 
formation about their decay, we assume that it looks like the decay of Z° . Theoretical and 
experimental arguments show that leptonic and hadronic decay channels of Z° have a branching 
ratio of ~ 1/3 — 2/3 [^6|. As hadronic channel is dominant, here we only consider this mode. 
It maximizes the flux of nucleons which at present are the dominant observable at ultra high 
energies. 

To mimic the softening of energy spectrum due to multiple decay level, we assume that decay is 



similar to hadronization of a pair of gluon jets. Experimental data |67| as well as MLLA (Mod- 
ified Leading Logarithm Approximation) |^8|, LPHD (Local Parton-Hadron Duality) ( [^] and 
references therein) and string hadronization model predict a softer spectrum with higher 
multiplicity for gluon jets than for quarks. 

We use PYTHIA program [^] for jet hadronization. This program, like many other available 
ones, can not properly simulate ultra high energy events, not only because we don't know the 
exact physics at VS^^GeV scale, but also because of programming limits. For this reason, we 
had to extrapolate simulation results for Ecu ^ lO^'^eF up to Ecu = lO'^^eV. Fig.|| shows 
as an example, proton and photon multiplicity in hadronization of a pair of gluon jets. At 
middle energies, the multiplicity per log(£') is roughly constant. The same behavior exists for 



other species. This is a known shortcoming of present fragmentation simulations [71| (See also 
Appendix 1) and makes spectrum harder at middle energies. As we would like to study the 
maximal flux of UHECRs and their effects on high energy backgrounds, this problem can not 
change our conclusions. 

In the simulation, all particles except e^,p^,z^, P and 7 decay. We neglect neutrinos mass and 
for simplicity we assume only one family of neutrinos i.e. fg- 

Contribution of the stable species in the total multiplicity and the total decay energy is sum- 
marized in Table |l|. For all species, more than 99% of the total energy belongs to the particles 
with energies higher than 10^*^6^ and lO^^ey respectively for two masses considered here. Ap- 
parently the mass of UHDM has little effect on the composition of remnants. However, one has 
to admit an uncertainty about this conclusion which is a direct consequence of the uncertain 
behavior of the multiplicity spectrum as mentioned above and in the Appendix 1. 



2.2 Interactions 

We have included roughly all relevant interactions between remnants (except u — i' and z/ — 
elastic scattering) to the simulation either analytically or by using the results of PYTHIA Monte 
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Table 1: Energy and multiplicity contribution in remnants of WIMPZILLA. 







lU el/ 




lU eV 


Part 




iViUlLi. /O 




Mnilti % 

iVlLliLl. /O 


e± 


6.7 X 2 


9.7 X 2 


6.7 X 2 


9.8 X 2 


P± 


11.8 X 2 


1.4 X 2 


11.9 X 2 


1.4 X 2 




18.4 X 2 


28.3 X 2 


18.1 X 2 


28.3 X 2 


7 


26.2 


21 


26.6 


21 



Carlo. Previous works |57] [58| |£^] |£2[ either don't consider the energy dissipation or take 
into account only the first order perturbative interactions (except for p — 7 where a fitting is 
used for N — ^ ^ N — tt cross-section) . An early version of the present simulation studied 
only energy dissipation of UHECRs using PYTHIA without considering low energy processes 
and interactions of secondary particles. It found a higher lifetime for UHECRs than present 
work. 

The main reason for considering only first order interactions is that it is usually assumed that 
interaction with CMB is dominated by the minimal process i.e — 7 — > — vr. However, the 
CMB spectrum even on its peak spans over a wide range of energies where radiation correction 
and hadronization become important. For instance, at Ecm = '^GeV the mean multiplicity 
is ~ 15 in place of 5 (after pion decay) in the minimal interaction. Moreover, in the galactic 
medium, the IR and visible radiations are comparable with CMB and play an important role in 
the energy dissipation of protons of £^ ~ 10^^ — IQi^^eV . In extragalactic medium, the number 
density of background high energy photons with [E > leV) is larger than visible and near 
IR. Another factor which accelerates energy dissipation of protons in the interaction with high 
energy photons is energy loss of leading proton in p — 7 interaction. It increases with energy 
(See Figj2|) and results a higher dissipation rate. 

PYTHIA can not simulate processes with invariant CM energies Ecm = \/^, Ecm < 2GeV to 
Ecm < 4:GeV (the lower limit depends on the interaction). Consequently, for smaller energies 
we have included only perturbative, first order interactions using analytical expressions. Table | 
summarizes processes which are included in the simulation, the energy range of analytical and/or 
Monte Carlo calculation of the cross-sections, and the cuts used for removing singularities at 
small energies or angles in the case of analytical calculation (results depend somehow on these 
cuts specially in the case of moderate energy and angular resolution of our program). For energy 
ranges that PYTHIA has been used, in general we use default value of various parameters of 
the program as defined in PYTHIA manual. A number of parameters have been changed for all 
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Table 2: PYTHIA Parameters 



Selected processes 


lept.-lept, lept.-had., unresolved 7, lowpt, diffract. 
& elastic had.-had., Radiative correction 


Status parameters 


2'^''-order , continuous pt cutoff. 


Pt cuts 


i^CM = ^, Pt = U 


Number of flavors 
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processes e.g. to make unstable particles like mesons decay etc. They are listed in Table |2|. For 
some interactions, the value of a few other parameters are changed. They are also listed in the 
Table ^ From now on s, t and u are Mandelstam variables. 

At very high energies, Eqm > W^'^eV for nucleon-nucleon and Ecm > W^^eV for other 
interactions, PYTHIA becomes very slow and the number of rejected events increases rapidly. 
For these energies we perform a linear extrapolation from lower energies as explained below. 



3 Evolution 



We assume that non-baryonic Dark Matter is totally composed of slowly decaying UH particles. 
From (|2|) and (Q) below it is evident that width and fraction of UHDM in DM are degenerate 
and in the evolution equations, decreasing contribution is equivalent to increasing lifetime. 
Boltzmann equation for space-time and energy-momentum distribution of a particle i is ||7^ (We 
use units with c = h = 1): 

p'^d,f'^'\x,p) - {Tii^p'^pP - e^Fj^p-)^ = -{A{x,p)+B{x,p))f^^{x,p)+C{x,p) + 

V{x,p)+S{x,p). (1) 



A{x,p) = Tiuii. (2) 
C(x,p) = T^^^^^J^J dPjf^'H^^P,)^^- (4) 
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X and p are coordinate and momentum 4- vectors; f^'^\x,p) is the distribution of species i] rrii, 
Ci and Tj, are its mass, electric charge and width = I/tj, t, is the lifetime; CTij is the total 
interaction cross-section of species i and species j at a fixed s; = ^g^^dpdn 

Lorantz invariant differential cross-section of production of i in the interaction of j and k; gi is 
the number of internal degrees of freedom (e.g. spin, color); dp = We treat interactions 
classically, i.e. we consider only two-body interactions and we neglect the interference between 
outgoing particles. It is a good approximation when the plasma is not degenerate. It is assumed 
that cross-sections include summation over internal degrees of freedom like spin; ^ is the 
differential multiplicity of species i in the decay of j; T^p is the connection; Fy an external 
electromagnetic field; and finally £{x,p) presents all other external sources. A{s) is a kinematic 
factor l7|: 

MPiiPj) = {(Pi-Pj)'^ - m^mj)^ = ^((s -mf - - ^mjrri^)^ . (6) 
The quantity Aa presents the probability of an interaction. 

In a homogeneous universe f{x,p) = f{t, \p\) and in (|^) the term corresponding to interaction 
with external electromagnetic field is zero. Therefore, to have a consistent formalism for evo- 
lution of distribution of all species, we don't include the synchrotron radiation of high energy 
electrons in a magnetic field. 

We only consider the evolution of stable particles and slowly decaying UHDM. The term (^) 
concerns only UHDM. In (|), the only non-zero term in the sum is the decay of UHDM. We 
assume that stable species don't have any interaction with UHDM and corresponding interac- 
tion integrals in (|3|) and (^) are zero. Due to the very large mass of UHDM, its momentum 
is negligible and we can assume that in comoving frame it is at rest. This permits to use its 
number density ridm which is more convenient for numerical calculation. 
In a homogeneous universe the metric in comoving frame is: 

ds'^ = dt^ - a^{t)6ijdx'dx^ . (7) 

and with respect to local Lorantz frame @) to (P) take the following form (in the following the 
species index indicates one of the stable species): 

^^-W^-^P^ = ^i-B{t,p)f^^\t,p)+C{t,p)+Vit,p)). (8) 
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E ^ dM^ 

^■KQiP^ dp 



^i'^iP) = A„ ■ (10) 



(11) 

dudm 3d _ _ 

~r ^dm -I- dm^dm- \^^) 

dt a 



In and (11) s depends on the angle between species j and k, and j and i. Consequently, it 
is not possible to use cross-sections integrated over angular variables. 
Evolution of a{t) is ruled by Einstein equation: 

d^ SvrG^ A , , 

^°(*> = Y.jS^ j ''PP''Ef">(t,p). (14) 

In a homogeneous cosmology, T*^*^ in Local Lorantz frame is the same as comoving frame and 

'~rii — „ — 2rnii 

^ comov " Loc.Lor.' 

Equations (^ and ([l^) determine the cosmological evolution of species. Due to interaction terms, 
even in a homogeneous universe, these equations are non-linear and coupled. It is not therefore 
possible to solve them analytically. Evolution equation for DM can be solved analytically for a 
short period of time. For other species if in (|8|) we consider absorption and production integrals 
as t and p dependent coefficients of a linear partial differential equation, (^) can be solved 
analytically ||7^. Giving the value of f^^\t,p), ndm{t) and a(t), at time t, we can then determine 
a{t + At), ndm{t + At), B{t,p), C{t,p) and V{t,p) for a short time interval At. /(*)(t + At,p) 
would be obtain from solution of partial differential equation (^) using difference method. The 
solution of metric and distributions in one step of numerical calculation are the followings: 

a(t + At) = a(t)exp(At(^Too + |)^). (15) 

ndm{t + At) = ndm{t) o.^x exp( ). (16) 

a'^(t) T 

f^'\t + At,p) = if^'\t,p') + At{C{t,p')+Vit,p')))eM-S{t,p')At) (17) 

p' = ^-^P^. (18) 
a(t) ^ ^ 

This prescription is more precise than a pure numerical calculation using e.g. difference method. 
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4 Numerical Simulation 



What makes numerical calculation of ( p^ ) to ( [ITD difficult is the extension on roughly 34 orders 
of magnitude of energy from W^^eV (radio background) to W'^'^eV (mass of UHDM) (from now 
on we call this energy range TZe)- Physical processes in this vast energy range have varieties of 
behavior, resonances, etc. Moreover, species have distributions which are orders of magnitude 
different from each others. In other term, these equations are very stiff. Semi-analytic method 
explained above helps to increase the precision of the numerical calculation. However, the 5- 
dimensional integration in (|ll|) is extremely time and memory consuming and it is impossible 
to determine it with the same precision. 

In the following we describe in detail the numerical calculation of (|l5| ) to (|l7|), as well as cos- 
mological model, initial conditions and backgrounds which have been used. 

4.1 Multiplicity and Cross-section 

For processes simulated by PYTHIA, we need to calculate total and differential cross-sections 
(see (^) and (0)). The former is given by the program itself. To determine the latter, we divide 
TZe to logarithmic bins (one per order of magnitude) and classify particles according to their 
momentum. The angular distribution of produced particles with respect to the axis of incoming 
particles in CM also is divided linearly to 90 bins. At a given s, the cross section in each bin 
Aaij = (JtotNij /N . Nij is the number of particles of a given species in the bin ij. N is the total 
number of simulated events. 

The same procedure is used for determination of In this case it is not necessary to consider 
the angular distribution because in the rest frame of WIMPZILLA the decay has a spherical 
symmetry. 

Because PYTHIA can not cover the totality of the energy range, for high energy bins we use 
a linear extrapolation in log p. The contribution of these energies i.e. EcM > W^GeV on the 
evolution of species is nevertheless small because the density of concerning particles is very low. 
The reason for adding them is not to have an artificial cut in the calculation. 
As mentioned above, for most processes including p— 7, PYTHIA can not simulate the interaction 
with Eqm < AGeV . p — 7 is the most important process for the energy dissipation of protons 
specially in this uncovered energy range where interaction with CMB photon is concentrated. 
In these energies we use directly the total cross-section obtained from experience |6^. For 
differential cross-section, we extrapolate angular distribution from higher energies and normalize 
it to the exact total cross-section. 
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4.2 Evolution Equations 



In (P) and (|ll|), the integrals over angular degrees of freedom can be separated from energy 
integrals and they don't depend on any cosmological or DM parameter. It is therefore very 
convenient to calculate them separately. 

We divide each 180° interval to 9 bins and use trapezoid method for integration. For the single 
integral in @ a better resolution with 90 bins has been used. However, our tests show that 
even the moderate resolution of 9 bins gives, up to a few percents, the same results as the more 
precise integration. This is a reassuring results and means that the triple integrals in production 
term also must be enough correct even with a low resolution. 

Calculation of production term for first-order interactions is more complicate. They are processes 
of type 2part. 2part. In the CM frame where all cross-sections are determined, analytically or 
numerically, the incoming and outgoing particles have the same momentum. This is equivalent 
to having a Delta function in the integrand. The numerical realization of this function specially 
with a moderate resolution is very difficult. Consequently, one has to analytically absorb this 
function into integrand. Details of the calculation can be found in Appendix 2. 
Numerical solution of evolution equation itself needs much better energy resolution due to stiff- 
ness of distributions. The resolution must be at least comparable to smallest quantities. Our 
tests show that a division to 680 logarithmic bins of the energy range TZe (i-e 20 bins per one 
order of magnitude) gives an acceptable compromise between precision and calculation time. 
We divide the interval z = [z^gc" 0.001] to 30 logarithmic bins and the last step is from z = 0.001 
to z = 0. The program is written in language and is highly modulable. It can be requested 
from the author. 

To test the precision of our numerical calculation, we have run the program without interaction 
terms. The error on total T^^ is 0.7% and on non-baryonic DM is practically zero. For other 
species it is 5.5% to 7.5%. Including interaction terms but not the decay of DM gives the same 
answer. This test is crucial for correct simulation of thermal equilibrium of the Universe. 



4.3 Cosmology Model and Initial Conditions 

We consider a flat universe with present value of parameters as the followings: f^M = 0.3, 

= 0.7, h = Ho/WOkmsec-^ Mpc'^ = 0.7 and Qb = 0.02/i-2. 
We fix the decoupling redshift at z^ec = 1100. The distribution of species at that time was 
thermal with a temperature T^ec = Tcmhi^dec + 1) = 0.26eV, Tcmh = 2.72SK [77] for e^, and 
7 and ^T^^ec for u and u. 
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In the same way, one can determine the temperature of Dark Matter Q at decouphng [78|: 



2 



dec 3 rp2 

Tdmizdec) = -j^ ^ (19) 

9*s^^ dm— dec 

m—dec is the decoupUng temperature of the Dark Matter. If we assume that Tdm— dec ~ l^^^eV, 
Tdm{zdec) < lO^^^eV. Therefore, our approximation Tdm = is quite justified. 
We know the density of species at present. Their initial value at decoupling depends on the 
equation of state of the Universe. However, it is exactly what we want to calculate! Consequently, 
we have to determine their value at decoupling approximately by neglecting the effect of UHDM 
decay. As the lifetime of WIMPZILLA is assumed to be much longer than the age of the Universe, 
the present value of densities after evolution must stay very close to our initial assumption. 
We define the initial densities as the followings: 

"-7 = ncoBE{zdec + 1) n,, = np = 3 X — X "-7 (20) 

np=^^^{zdec + '^f n^-=np np = = Q (21) 

Tflp 

p^{VlM-^hot-^b{l + ^)) „ 7r2 , 

ndm = ^{Zdec + lf nnot = -g.T^ml,- (22) 

nidm 30 

Knowing initial density and temperature, other quantities like chemical potential and distribu- 
tions can be determined ||7^. We assume that the age of the Universe tq = 14.8Gyr. This 
quantity also depends on the equation of state and we use it only for fixing the lifetime of 
WIMPZILLA. 

4.4 Backgrounds 

Apart from CMB and relic neutrinos which are included in the initial conditions, we don't in- 
clude any other background to high redshift distributions. For z < 3, we add near-IR to UV 
emissivity of stars [79| |8C] to equation (^). Far-IR which is very important for energy dissipa- 



tion of UHECRs, is not added because there is very little information about its evolution with 
redshift. High energy backgrounds are not added because we want to be able to distinguish the 
contribution of remnants from other sources. 

Adding backgrounds by "hand" evidently violates the energy conservation of the model, but 
there is not any other simple alternative method. Moreover, the violation is very small and 
comparable to numerical errors. 



"^One should consider Tdm as an estimation of kinetic energy scale rather than a real temperature because if 
ultra heavy particles exist, they could never be thermalized. 
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Here a comment is in order: why have we included star backgrounds and not the synchrotron 
radiation ? Star background is considered as an external source. By contrast, synchrotron radi- 
ation concerns high energy electrons which are involved in the evolution. In (Q), the elimination 
of interaction with an external field in a homogeneous universe means that the probability for 
production and absorption of synchrotron photons is the same. It is not therefore possible to 
consider their production without their absorption i.e. interaction of electrons with external 
magnetic field. 



5 Results 

Figure |^ shows the energy flux of high energy protons and photons in a homogeneous universe. 
Fig.^ shows the same quantity for all species. The GZK cutoff is very transparent. With our 
background model it begins at ~ lO^^^'^eV for protons due to p — 7 interaction (see optical 
depth in Fig.^) and at ~ lO^'^eF for photons due to production. For purely kinematic 
reasons, proton cutoff is much shallower than photon one. The resonance of p — 7 interaction is 
very close to p rest mass where A{s) in (^) is very small. Moreover, division of by E reduces 
the effect of absorption on the distribution. 

According to Fig.^, even the shortest lifetime we have considered, can not explain the observed 
flux of protons (in contrast to |^2| that assumes a 2-particle decay mode). The same figure 
shows also the flux without energy dissipation. It is compatible with Ref. ||6^ which assumes a 
hadronic decay but does not consider the energy dissipation of secondary particles. In this latter 
case, the lifetime must be ~ 4—6 orders of magnitude larger than the age of the Universe. On the 
one hand, this result proves the role of a realistic model of energy dissipation in the estimation 
of mass and lifetime of UHDM. On the other hand, it shows the importance of clumping of Dark 
Matter, i.e. most of observed UHECRs must come from nearby sources. This conclusion is 



independent of the source of UHECRs. A similar conclusion has been obtained in |58| by fitting 
the expected flux from extragalactic sources and from the Halo on the data. However they don't 
consider the dissipation. In the next section we show that even in a halo, the dissipation of 
photons energy is significant. 

In the case of decaying UHDM hypothesis, the most important source is the Galactic Halo. 
Before trying to make a simple model of halo in the next section, we discuss some of other 
conclusions one can make from this study. 

From Fig.^ and by taking into account the fact that p — 7 cross-section is ~ 10^ times smaller 
than p — p, with present statistics of UHECRs, less than one photon shower could be observed. 
The comparison of EGRET data for 10*^6^ < E < lO^^eV with our calculation shows that it 
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is compatible with a short hfetime UHDM. Future observations of GLAST at E > W^^eV is 
crucial for understanding the source of UHECRs because one expects an additional contribution 



from synchrotron radiation of ultra high energy electrons at -E ~ 10^^ — lO^'^e^ [^] |62]. 
Fig.^ illustrates the total and partial optical depth, i.e. 13{t,p) in equation (|9|) at z = 0. It 
shows that for protons with E ~ W^^eV, even a source at ~ 20Mpc must be very strong to 
be able to provide the observed flux because its flux is reduced by 11 orders of magnitude! 
The optical depth of protons must be even larger than what we have obtained here because we 
didn't take into account the far-IR and radio backgrounds. This makes difficulties for the recent 



suggestion by Ahn E.J., et al. [40| that Virgo cluster can be the only source of UHECRs, even 
if the magnetic field of Galactic wind is as strong as what is considered in that work. The main 
challenge is finding a conventional source with enough emissivity of CRs at ultra high energies. 



5.1 CMB Distortion and Entropy Excess 

To see if a relatively short living UHDM can distort CMB, we reduced the spectrum of photons 
for a decaying DM from the spectrum with a stable DM. FigJ^ shows the result for z = Q. There 
is no distortion at least up to 1 to 10^ parts for E < ?>eV . However, one should note that this 
conclusion depends somehow on the cross-section cuts at low energies. Nevertheless, giving the 
fact that CMB flux at its maximum is ~ lO^'^ times larger than the rest of the spectrum, it 
seems very unlikely that it can be disturbed, otherwise SZ effect due to stars had to be observed 
everywhere. The distortion of CMB anisotropy will be studied elsewhere. 

We have also examined the entropy excess separately for each species. There is no entropy 
enhancement except for e"*" and p~ which in our model are absent from the initial conditions. 
Comparing to other species, their contribution is very small and negligible. 



5.2 Baryon and Lepton Asymmetry Generation 



It has been suggested |81] [22] that the decay of UHDM may be able to generate additional 
baryon and lepton asymmetry. At GUT scale, i.e. the mass scale of WIMPZILLA, we expect 
such processes and the fact that a late time decay is out of thermal equilibrium and satisfies 



Sakharov conditions for baryogenesis [78| make growing asymmetry plausible. 

The rate of baryonic (or leptonic) number production by decay of UHDM in comoving frame 

can be expressed as: 

^^^V^ + - n,) = (23) 

at a{t) T 
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e is the total baryon number violation per decay. The solution of this equation is: 

AK-ng) = enrf„(to)(l-exp(-^))i^±^. (24) 
. ^(^^ = |!^(l_exp(-^))il±^ (25) 

If to = tdec, i^p^ ~ 10~^^ (for rudm = W'^'^eV). Therefore AB ~ IQ-'^^e at z = 0. As e can 
not be larger than total multiplicity, ~ 1000, AB < 10^^^, i.e. much smaller than primordial 
value ~ 10-1°. 

We tested this argument by assuming e = 0.1 at all energies, i.e. etot = O.lAitot- Evidently np is 
smaller, but the change of Up is too small to be measured. The same is true for number density 
of leptons, but energy density of leptons with respect to anti-leptons increases by an amount 
comparable to e. 



5.3 Equation of State of the Universe 

Decay of UHDM gradually changes part of CDM to HDM and thus changes the equation of state 
of the Universe. Fig.^ shows the variation of equation of state. For r ~ SOtq or larger, it would 
be too small to be measurable. For smaller r, DM decay plays the role of a running cosmological 
constant. A complete study of this issue and comparison with data is under preparation. 



6 Halo 

To see the effect of clumping of a decaying UHDM on the flux of UHECRs, here we try to make 
a very simple model. A complete treatment of halos will be reported elsewhere. 
We consider a halo as a uniform over-density with a limited size at z = 0. This is simulated 
by following the decay of UHDM and evolution of remnants for a time comparable to the 
propagation time in the halo. Evolution equation and energy binning is taken to be the same as 
in the homogeneous universe case with a{t) = cte. Because we want to study the propagation 
of remnants in a volume comparable to the Galactic Halo, we consider time steps equivalent 
to Wkpc. Our tests show that after a few steps (~ 7), the accumulation rate of ultra high 
energy particles becomes very slow. We consider two cases. In the first case we simply evolve 
distributions for a number of steps (up to 30). In the second case, after some steps, we stop the 
decay of the Dark Matter to simulate an inner halo of MACHOs. Then, the evolution is continued 
for more 5 steps (i.e. 50kpc) to simulate propagation through MACHOs. Evidently this 
model is very approximative. We use it only to make a crude estimate of production and 
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absorption of UHECRs. 



6.1 Initial Conditions and Galactic Backgrounds 

Galactic baryonic matter is taken to have a thermal distribution with = W^K. We assmne 
that baryonic over-density is biased with respect to DM i.e. the fraction of baryons to DM is 
larger than its mean value in the Universe. With these assumptions the initial number densities 
can be expressed as: 

bSpc 

— Tip— — 



mp 



Up = = 
(1 - b)Spc 

ridm = (26) 

rridm 

b is the fraction of baryons in the halo. In the following b = 0.3, i.e. ~ 2 times primordial 
value [p5| in the cosmological model explained above. 6 is the mean over-density of the Halo. 



Inspired by universal halo density distribution of NFW [86], we consider a halo with characteristic 
radius (i.e. virial radius) r2oo = 0.12Mpc. According to NFW distribution and by definition 
this means 6 = 200. These parameters defines a halo of mass Mh = 6 x 

Neutrino density is assumed to be the same as relic at z = 0. For photons, in addition to CMB, 
we consider a galactic background as the following: 

Galactic IR and visible backgrounds are not very well known. We use the results of the model 
developed by DIRBE group for detection of extragalactic component of the IRB ||8^. We 
consider the observed value of IRB after elimination of Inter-Planetary Dust (IDP) contribution 
as the galactic background. It is just an estimation of average galactic IRB. It is not clear if we 
can extend the local value of IRB to whole galaxy or take it as a representative average. For 
this reason we also increase it 10 times (probably an extremely high value) to see the effect on 
the energy dissipation of UHECRs (see below for conclusions) . Our simulation does not include 
radio background. 

For soft and hard X-Ray galactic backgrounds, we use the model developed for extraction 
of extragalactic component from ROSAT and ASCA observations |^8|. It considers GXB as 
two thermal components, a soft component with Tg^ = 70eV from Local Bubble, and a hard 
component with Thx = 14:5eV from hot gas, probably in the Halo. We add also the extragalactic 
component for 0.25keV < E < lOkeV. 
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7 Results 



Fig.^ shows the distribution of high energy protons and photons for a uniform halo and for a 
hafo that its inner part is composed of MACHOs. Only the result for mdm = W^'^eV with 
r = 5tq and r = SOtq is shown. Because of importance of the spectrum close to observed trough 
in the interpretation of the results, Table Q summarizes the numerical value of simulated and 
observed spectrum. 

For photons the trough of GZK cutoff is shallower than in a homogeneous universe and for 
protons it is practically absent. Consequently, the calculated flux at ~ W^^'^eV is somehow 
higher than observation. However, at higher energies simulation results specially one with a 
MACHO halo is in the 1 — a error range of the observations. 

The simplest explanation of having a smaller observed flux close to the minimum of the spectrum 
can be the need for increasing the lifetime somehow (but not by many orders of magnitude as 



suggested by previous works |60] |6^). However, the lack of a minimum in the simulated 
spectrum]^ means that in some way our simulation is not exact. It can be due to a too simple 
halo model. Other possibilities are a different (probably harder) decay spectrum and/or energy 
loss in the Galactic magnetic field. If these suggestions are true, fluxes will be smaller and there 
is no need for increasing the lifetime. A more complete simulation of the Halo and magnetic 
fleld and a better understanding of non-baryonic DM distribution is necessary for making any 
definitive conclusion. 

The minimum in the spectrum can also be interpreted as a wider distribution of sources. In this 
case, it is hardly probable that sources responsible for CRs at lower energies can explain this 
behavior because even if the spectrum of these sources can be extrapolated with the same slop 
to higher energies, it can not explain the rising slop of the spectrum. 

As the IR background is crucial for energy dissipation of UHECRs, we have also increased it to 
10 times of the DIRBE model to see the effect. Proton fiux at high energies slightly decreases, 
but it can not explain the observed minimum. We have also tested the distortion of the CMB 
by remnants as described for a homogeneous universe. There is no distortion up to at least 1 to 
10*^ for E < WeV. 

Summarizing this section, it seems that a decaying UHDM with a lifetime as short as r ~ 
10 — IOOto is not ruled out by present observations. 



'in fact optical depth has a maximum at E ~ 3.3 x lO^^eV. 
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8 Conclusions 



The main purpose of this work was showing the importance of a reahstic physical model for 
finding the answer to the mystery of UHECRs and introducing readers to the program that has 
been developed for achieving this goal. 

We showed that TZ = ^tq/t ~ 0.1 — 0.01 where ^ is the contribution of UHDM in the Dark 
Matter. This value is larger (or equivalently the lifetime is shorter for the same contribution) 
than what has been suggested in previous works which don't consider the dissipation |^8| []6C| ]. 
If a more realistic halo model confirms this conclusion, the cosmological implication of a UHDM 
can be important. 

We have studied a special decay mode. Any other mode that produces invisible WIMPs or 
leptonic or semi-leptonic modes decreases the lifetime of UHDM. If more nucleon are produced, 
the lifetime must be longer but it seems less probable than other cases. 

We showed the reciprocal influence of UHECRs and backgrounds on each others. Consequently, 
whatever the source of UHECRs, it is extremely important to correlate their observations to the 
observation of high energy photon and neutrino backgrounds. 

This work is the first step to a comprehensive study of the effects of a decaying Dark Matter. 
Other issues like a realistic model for halos, effects on the determination of cosmological param- 
eters and equation of state and comparison with more data from cosmic rays and high energy 
backgrounds remain for future works. 

Note: Shortly after completion of this work the Lake Baikal Experiment collaboration pO[ | have 
published their upper limit on the flux of high energy neutrinos. It is well above what is obtained 
in our simulation (see Fig.^). 



Appendix 1: Fragmentation in MLLA 

The MLLA treats fragmentation as a Markov process. Consequently, the differential multiplicity 
is proportional to splitting function P(x)(See e.g. [|6^] ): 

For a gluon jet and at xe ^ 1, P{xe) ~ and one expects that: 

dM{xE) 



d{ln{xE)) 



oc Total Num. of sphtting oc Ej. (28) 
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In PYTHIA however, at high energies E > lOOGeV, increasing CM energy just increases the 
probabihty of having more high energy fragments which simply escape fragmentation. A conse- 
quence of this behavior is the narrowing of the multiphcity distribution (See Fig.^) in contrast to 



theoretical prediction i.e. KNO |89] scaling for Ej ^ oo. This limit is obtained at parton level. 
However, if LPHD is valid, one expects the same type of behavior at hadron level. Another 
factor that can explain, at least partially, the deviation from theory is the decay of hadrons in 
our simulation. It increases total multiplicity more than its statistical variation and makes the 
distribution narrower, but it can not explain the absence of a tail of large multiplicity events. 
Even after decay, these events should keep their difference with average. As both MLLA and 



Monte Carlos fail to reproduce observations in relatively low energies [|63] [71], the exact behavior 
at high energies is not clear. 



Appendix 2: Production Integral for 2^2 Processes 

Cross-sections are Lorantz invariant. However, to have a unique expression for what are simu- 
lated and what are calculated analytically, we determine all of them in their CM. 



The triple integral in (11) depends on three momentum variables corresponding to the mo- 
mentum of two incoming particles and one of the outgoing particles that its evolution is under 
calculation. In the case of a 2 — > 2 process, the amplitude of the momentum of final particles 
depends only on s. If: 

p = {E,pcos(j)sin6,ps'm4)sm6,pcos6). (29) 
is the 4-momentum of outgoing particle, and p' is its counterpart in the CM: 

P' = ^ ^ • (30) 

p' = rp (31) 

m and m' are the mass of out-going particles, T is the boost matrix. The equality of two ex- 
pressions for p' leads to an equation that can be solved for one of the angular variables in ([29|). 
The calculation is tedious but strait forward. With respect to 0, the equation is 4*'* order but 
in the case of a homogeneous cosmology where the boost matrix depends only on the relative 
angle between incoming particles, it depends only on tan^(^) and is analytically solvable. In 
this way the integration over 4> in ( [III) reduces to sum of integrand evaluated at the roots of the 
equation. 

This method provides a general way to deal with this problem and is more convenient than 
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doing calculation for each cross-section separately. 
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Figure 1: Proton and photon multiplicity in hadronization of a pair of gluon jets for Ecm = 
IqU _ io24gy_ Yoi Ecm > W^^eV the curves correspond to extrapolation from lower energies. 
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Figure 2: Energy fraction of leading proton (dashed line) and mean energy of protons 
(solid line) with respect to the CM energy in p — j interaction. 
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Table 3: Interactions 



Interaction 


Analyt. Cal. 
EcM Rang 


Cuts 


PYTHIA EcM 
Rang 


PYTHIA Parameters 


Ref. Analytic 
Cross- Sec. 


— elastic 


All 


> 10-3, 
■!^| > 10-M-!i| > 10-2 




- 


El 


— e~ elastic 


£cA/ < 4GeV 




included in e"*" — e- 
— * . . . 


- 


11 
75 


e+ - e" -> 27 










11 


— ^ . . . 




- 


4GeV^ <EcAi< 10^ 


Init. Braamst. 


- 


p± - e± ^ . . . 
(all combinations) 












p±-p±^... 
(all combinations) 






3GeV < EcM < 10^ 






V — e~ elastic 


All 








Izl 


V — elastic 












V — ^ . . . 






4Ge\/ < EcM < 10^ 


- 




V — . . . 












V — ^ . . . 






4GeV <EcM < Iff 






u — p^ . . . 


- 


- 


11 


- 


- 


u — V 


- 


- 


11 


- 




7 — e± ^ 7 — e± 


-Be A/ < 4GeV 


> lo-^ 

III > 10-2, Ifl > 10-2 






J L 


7 — e± . . . 






4GeV < EcM < 10« 






7 - p± ^ . . . 


lGe\/ < £ < 
4GeF 






Ecm""" = IGel^, 
Pt™"" = 0.5GeV, 
Singul. Cut = 0.25GeV 


il 


7 — 7 — 


EcM < 3GeV 


III > 10-2, Ifl > 10-2 






ii 


7 - 7 ^ . . . 






3GeV <EcM < 10^ 
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Table 4: Energy flux of UHECRs close to GZK cutoff. Note that there is not a complete 
agreement between estimated flux by different Air Shower detectors. The values presented here 
are mostly based on AGASA data before 1998. 



log(£;)eV" 


r = 5ro, logE^J{E) 
sec""' sr~^ 


r = 50to, \ogE^J{E) 


Observed log E^J{E) 

m^^ sec"""^ sr~^ 


18.525 


5.4 


4.4 


6.15 ±0.01 


19.025 


5.9 


4.9 


5.5 ±0.05 


19.525 


6.35 


5.3 


5.08 ±0.2 


20.1 


7 


6 


5.4 ± 1.4 


20.525 


7.36 


6.3 


~ 5±2 


21.025 


8.09 


6.8 






Figure 3: Energy flux for protons and photons. Solid line mdm = W eV, r = Btq, dot line is the 
spectrum without energy dissipation for the same mass and lifetime, dashed line mdm = W^'^eV, 



T = 50to, dash dot mdm = 10 eV, r = 5tq, dash dot dot dot rridm = W eV, r = SOtq. For 
protons, data from Air Showers detectors |13| is shown. Data for photons are EGRET whole 
sky background |82] and upper limit from CASA-MIA |83|. 
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Figure 4: Energy flux of stable species. Description of curves is the same as Fig. 
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Figure 5: Total optical depth of species and contribution of backgrounds. Solid line is total 
optical depth, dot line contribution of 6=*=, dashed line p"^, dash dot u&Liy, and dash dot dot dot 
7. Dependence on lifetime and mass of UHDM is negligible. 



2,0x10 - 




-10 -8 



Figure 6: Fraction of distortion in photon distribution with respect to a stable DM. 
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Figure 7: Energy density of the Universe. Solid line a stable DM, dash line r = 50ro, dash dot 
r = 5ro. Dependence on the mass of UHDM is negligible. 
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Figure 8: Flux of high energy protons and photons in a uniform clump. TfidYn — 

W^eV, T = 5ro 

and T = SOtq. Dash dot and dash dot dot dot lines presents UHDM halo. Solid and dashed 
lines show a halo of UHDM and MACHOs. Data is the same as in Fig.^. For protons the effect 
of increasing lifetime of UHDM is more important than presence of MACHOs. Photons trough 
is more sensitive to presence of MACHOs. 
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Figure 9: Distribution of multiplicity in hadronization of a pair of gluon jets. Solid line is the 
highest energy Ecm = W'^^eV; Long dash is lowest energy EcM = 10"ey. KNO distribution 
is shown also (solid line). 
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